Onsager approach to ID solidification problem and its relation to phase field 

description. 
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We give a general phenomenological description of tlie steady state ID front propagation problem 
in two cases: the solidification of a pure material and the isothermal solidification of two component 
dilute alloys. The solidification of a pure material is controlled by the heat transport in the bulk 
and the interface kinetics. The isothermal solidification of two component alloys is controlled by the 
diffusion in the bulk and the interface kinetics. We find that the condition of positive-definiteness 
of the symmetric Onsager matrix of interface kinetic coefficients still allows an arbitrary sign of 
the slope of the velocity-concentration line near the solidus in the alloy problem or of the velocity- 
temperature line in the case of solidification of a pure material. This result offers a very simple 
and elegant way to describe the interesting phenomenon of a possible non-single-value behavior of 
velocity versus concentration which has previously been discussed by different approaches. We also 
discuss the relation of this Onsager approach to the thin interface limit of the phase field description. 
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Introduction. In the recent years the phase field ap- 
proach to solidification problems has attracted the at- 
tention of many researches (see, for example, [l[ and ref- 
erences therein) . It was originally introduced as a math- 
ematical tool to solve the free boundary problem without 
directly tracking the interface position. However, more 
recently it has also been considered as a physical model 
which can bring additional information compared to the 
sharp interface approach. In particular, it was observed 
that the general believe, that steady state ID front prop- 
agation with positive velocity 



V = Vo{At - 1) 



(1) 



is possible only if {At — 1) > (see, for example (H), 
is not the general situation. Here V, is the steady state 
front velocity, Vq is the characteristic velocity which is 
proportional to the kinetic growth coefficient; At = 
(Tjvf — Tq)cp/L is the dimensionless undercooling, Tm 
is the melting temperature and Tq is the temperature in 
the original phase far away from the interface; Cp is the 
heat capacity which is assumed to be the same in both 
phases; L is the latent heat. Karma and Rappel (KR) 
[3| introduced the thin interface limit of the phase field 
description and found that 



V = 



Vo{At - 1) 



(2) 



where Dt is the thermo-diffusion coefficient, a is a posi- 
tive number of order unity which depends on the details 
of the model, and W is the thickness of the interface in 
the phase field description. In the phased field model 
discussed in Q there is no any restriction on the pa- 
rameter VqW/Dt and the velocity may be positive for 
{At — 1) < . The same result was obtained for the 
isothermal solidification of alloys by many authors start- 
ing from a paper by Lowen et al. |j| in the framework of 
phase field description and also by Aziz and Boettinger 
who use a more phenomenological approach. In the 
case of alloys the deviation from equilibrium is defined 



by Ac = {Cl - Co)/{Cl - Cs) instead of At. In the 
two phase region of the phase diagram < Ac < 1. 
Here Cl and Cs are the equilibrium concentrations of 
the initial and growing phase, respectively, and Co is the 
concentration of the initial phase far from the interface. 
They found that the steady state growth is possible also 
inside of the two-phase region of the equilibrium phase 
diagram. 

From the numerous papers on the derivation of the 
sharp and thin interface limits from a phase field model 
we should also mentioned the work by Elder et al. 
0, Umantsev Q and also the paper by Korzhenevskii, 
Bausch and Schmitz Q which contain many details and 
technical points. The basic results of all these descrip- 
tions have the structure of Eq. ([5]) in the vicinity of 
{At{c) — 1) <C 1 and eventually lead to the non-single- 
value behavior of the velocity as a function of the driving 
force in the case of a negative "kinetic coefficient", Fig.l. 
In this case the branch which is described by Eq. ([2]) 
(dotted line in Fig.l) is linearly unstable (sec, for exam- 
ple. Hi, 0, II) while the "high velocity" branch of the 
mentioned non-single-value behavior is linearly stable. 

Qualitatively the same results have been obtained by 
the numerical solution of ID motion of the atomically 
rough interface in binary alloys . In this model instead 
of the phase field order parameter the authors used the 
fraction of the atomic places which belongs to the grow- 
ing phase. This fraction changes from to 1 during the 
growth. The evolution equations for this quantity to- 
gether with the concentration fields in the two phases 
are given by Eqs. (5.1)-(5.3) in The numerical anal- 
ysis of has shown that both types of curves in Fig.l 
are possible. However, the unstable (dotted line) branch 
was not seen in this dynamical simulations. 

The purpose of this Communication is to give a gen- 
eral phenomenological description of the steady state ID 
front propagation problem in two cases: i) the solidifica- 
tion of a pure material which is controlled by the heat 
transport in the bulk and the interface kinetics; ii) the 
isothermal solidification of two component dilute alloys 
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FIG. 1: Schematic dependence of the steady state velocity V 
vs. the dimensionless undercooling At. The curve 1 corre- 
sponds to the case aWVo/Dr < 1 while the curve 2 corre- 
sponds to the opposite case, aWVo/ Dt > 1. 



which is controlled by the diffusion in the bulk and the 
interface kinetics. Describing the interface boundary con- 
ditions we use only the general phenomenology of linear 
non-equilibrium thermodynamics in the spirit of the On- 
sager matrix of kinetic coefficients which has the proper 
symmetry and is positive-definite as required by the sec- 
ond law of thermodynamics. This approach does not 
assume any specific model of the interface and makes no 
assumption on its thickness. The only requirement, as 
for any macroscopic theory, is that the thickness is small 
compared the macroscopic lengths. We will see that two 
mentioned problems arc very close to each other and can 
be formally mapped onto each other. The mentioned 
restrictions on the Onsager matrix of kinetic coefficients 
are not sufficient to determine the sign of the slope of the 
velocity-concentration line near the solidus in the alloy 
problem (or of the velocity-temperature line in the case 
of solidification of a pure material). This result offers a 
simple way to describe the mentioned above phenomenon 
of a non-single- value behavior, Fig.l. 

The sharp {W — > 0) and the thin interface limits of 
the phase field description should lead to the effective 
macroscopic description with the boundary conditions in 
the spirit of Onsager relations, where the elements of 
the Onsager matrix arc expressed in terms of the phase 
field model parameters. Indeed, these limits really cor- 
respond to such a description. However, the mentioned 
condition of positive-definiteness of the matrix of kinetic 
coefficients turns out to be a nontrivial issue for the thin 
interface limit and will be discussed in more details. 

Growth of a pure material with heat transport. We 
assume that phase 1 grows at the expense of phase 2 by 
a ID front propagation with the steady state velocity V. 
In the bulk we have the thermal-conductivity equation. 



In order to write down the boundary conditions at the 
interface we follow the description and notations given in 



{p-2 - iJ-l)/TM 



M 



= AV 
= BV- 



-BJe 
CJe , 



(3) 
(4) 



where is the chemical potential of the corresponding 
phase i at the interface. According to the energy conser- 
vation at the interface, flux Je is defined by Eqs. (51)- 
(52) in 
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AiVTi = VTmSi^Je, 
-A2VT2 = VTmS2~Je 



(5) 
(6) 



Here 5*1(21) and S'2(T2) are the entropies of two phases 
and Xi is the thermoconductivity of the phase i. The 
elements of the Onsager matrix, which is symmetric 
and positive-definite, obey the conditions A,C > and 
B^ < AC. Rk = CTli is the Kapitza resistance and the 
cross coefRcient B describes the way the two entropies 
are shared between the two sides of the interface during 
growth (for a more detailed discussion of the physical 
meaning of the different Onsager coefficients in this case 
see [HI). 

For the steady state one-dimensional problem VTi = 
and Ti = To -t- L/cp where L = Tm[S2{Tm) - Si{Tm)] 
is the latent heat and Cp is the heat capacity; Tq is the 
temperature in the original phase far away from the in- 
terface. We note that in order to obtain the relation 
Ti ^Ta + L/cp one should expand the entropies near the 
equilibrium temperature Tm in the energy conservation 
condition (AiVTi - AzVTs) = VTm[S2{T2) - S'i(Ti)]. 
Now expanding the difference of the chemical potentials 
near the equilibrium temperature Tm, we find 

M2(r2)-/il(ri) = {S2-Si){Tm-Ti) + S2{Ti-T2) , (7) 



and finally we get 



V 



-CTISiS2+BTm{Si+S2)] 



(8) 



where At = {Tm — To)cp/L. We have used the usual 
notation for solidification of pure materials. We see that 
the sign of {At — 1) in general is not determined by the 
Onsager restriction B'^ < AC. However, it is well defined 
in two cases: i)B = and ii) in the "isothermal" case, 
Ti = T2. In the later case the growth rate is controlled 
by the "isothermal" kinetic coefficient which is strictly 
positive due to the mentioned restriction, B^ < AC 11| : 



V 



(M2 - pti) 



L^{At - 1) 



TmA[{1 - By {AC)] 



By{AC)] 



(9) 

Karma and Rappcl obtained in their thin interface 
limit B = C ~ and a coefficient A which may even 
be negative (/? in their notation). They discussed this 
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"counterintuitive" issue and gave some natural explana- 
tion for this phenomenon. We will return to this point 
later. 

Isothermal alloy solidification in the dilute limit. We 
discuss the steady state propagation of a ID front with 
velocity V during solidification of a two component alloy 
at a given temperature T. The concentration of B atoms 
is Ci{x) in phase 1 and C2{x) in phase 2. In the bulk 
these concentrations are described by diffusion equations 
with diffusion coefficients Di and D2- In order to write 
down the boundary conditions in this case we use the 
same phenomenological approach but adapted to the al- 
loy situation. Onsager relations connect two fluxes Ja 
and Jb (at the boundary) of atoms A and B to two driv- 
ing forces SfiA and dfiB which are usual differences in 
chemical potentials at the boundary. While the bulk is 
described by diffusional equations for the concentration 
fields for each phase, we still need three boundary con- 
ditions at the interface. One is the conservation of B 
atoms at the interface. We have also to relate the two 
concentrations Ci and C2 on both sides of the interface 
to the growth velocity and gradients of the concentra- 
tions. In the equilibrium these two concentrations are 
just the liquidus and solidus concentrations. When the 
velocity is finite these two concentrations deviate from 
the equilibrium values. We write (see, for example, [l2| 
and references therein): 



SfiA/T = AJa + BJb , 
SfiB/T = BJa+CJb , 



(10) 
(11) 



This Onsager matrix should be positive-definite: A and 
C must be positive and < AC. According to the con- 
servation of B atoms at the interface we also have [iGl 



L»iVCi = VCi - Jb 

-Z?2VC2 = VC2-JB 

V = Ja + Jb- 



(12) 
(13) 
(14) 



Eq. (14) is written for substitutional alloys. For inter- 
stitial alloys V = Ja- For dilute alloys the chemical 
potential are [Tsj 



S^ia/T 
S^ib/T 



{Ci ~~ C2) - 
HC2/C1) 



(Cl - Cs) 
-HCs/Cl) 



(15) 
(16) 



Here phase 1 grow at the expense of phase 2. Ci and C2 
are the concentrations of B atoms at the interface and Cs 
and Cl are their equilibrium values; {Cl ~Cs) ^ {Tj^j — 
T) /T is proportional to the deviation of the temperature 
from its equilibrium value for a pure A material. Di and 
D2 arc the diffusion coefficients. 

According to the mass conservation at the interface for 
the steady state ID growth wc have Ja = V{1 — Ci) and 
Jb = VCi because there is no gradient in the growing 
phase 1. These relations are written for the substitutional 
alloys. For the interstitial alloys Ja = V. However, in the 
dilute limit there is no difference between these two alloys 
because Ci ^ 1 and can be neglected in the expression 



for Ja for the substitutional alloys. Moreover, the global 
mass conservation requires that Ci = Co, where Co is 
the concentration in original phase 2 far away from the 
interface. Solving the resulting system of equation we 
find the transcendental relation between velocity V and 
the initial concentration Co: 



In 



Cs_ 

Cr 



Cl-Cs 
Co 



-V{A/Co + B) =V[B + CCo]. 

(17) 

If the concentration Co is close to Cs and the velocity V 
is small we find, expanding logarithm up to linear order 
in (Co - Cs) and V, 



V 



[Cl - Cs){Cs - Co) 
Cs[A + CClCs + B{Cl + Cs)] 



(18) 



For the general case of not dilute alloys this equation 
reads 



V 



[fi'{Cs)/T]{CL~Cs)iCs~Co) 



A{1 - Cl)(1 - Cs) + CClCs + B[{Cl + Cs) - 2ClCs] 

(19) 

where /"(C) is the second derivative of the free energy 
/i(C) of the growing phase 1 with respect to the con- 
centration. From this expression it is clear that in the 
presence of the cross coefficient B the sign of (Cs — Co) 
is not determined by the condition B^ < AC and also 
depends on C^ and Cs- Moreover, if the sign in the 
square brackets of Eq. is negative and Co > Cs for 
small positive velocity V then wc find for Co = Cs apart 
from the solution V = second solution with positive V. 
If the expression in the square brackets is negative but 
small, we can expand the logarithm up to linear order in 
(Co — Cs) and up to quadratic order in V and find: 



(Cl - Cs)iCs - Co) 
Cs 



V[A + CClCs + B{Cl + Cs)] 
+V^[A + BCs]^/i2CL) (20) 



This expression shows that with increasing V the curve 
V = V{Co) first goes into the two-phase region, then 
turns back having another solution with finite velocity at 
Co — Cs and then goes into the one phase region (see 
Fig. I). Eventually, for Co — >■ the velocity, according to 
Eq. (HZl), becomes V ^ {Cl - Cs)/A (Tm - T)/{TA) 
as for the solidification of a pure material. 

First of all, we would like to mention the clear analogy 
between two discussed problems. From the basic equa- 
tions we see this analogy if we relate V Ja, Je ^ Jb, 
TmSi{2) ->■ (^1(2), S^i S^iA, ST/TIj Sub/T and 
apart from some thermodynamical prefactors { At — 1) — )■ 
(Cs — Co). Then, the case 5T = in the pure mate- 
rial problem corresponds to zero values oi 5^b- This, in 
turn, corresponds to a frequently used assumption that 
the partition coefficient k = C1/C2 equals to its equilib- 
rium value fco = Cs / Cl ■ In this case as in the pure ma- 
terial problem stationary growth is possible only in the 
one phase region of the phase diagram. Actually it seems 
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that this result is in agreement with the phenomenologi- 
cal description of [1] and also 

Discussion and conclusion: Thin interface limit of 
phase field models vs. Onsager approach. We discuss 
the thin interface limit using the KR description for the 
temperature field for a flat interface. The correspond- 
ing problem for alloys leads to basically the same results 
(see, for example, Ji., 8,, JJ,]). Originally it was designed 
to increase computational power of the method by us- 
ing larger values of the interface width W and to mimic 
local equilibrium boundary conditions Q. Let us have 
a closer look at this limit from more physical prospec- 
tives. In the thin interface limit of Q the tempera- 
ture distribution r(a;) close to the interface is given by 
Ti{x) = T{0) + GiX where Gi is the temperature gra- 
dient in the i-th phase {i — 1,2) at a; = 0. At a; = 
the temperature Ti = T2 = T(0) and in this description 
the Kapitza jump is absent, Ti — T2 = 0. One should 
note that the value of T(0) in this linear extrapolation 
procedure is different from the real value of the smooth 
temperature field at the middle of the interface obtained 
by the phase field simulations. The given linear ex- 
trapolation of the temperature field reasonably coincides 
with direct phase field results only for ^ |a;| ^ 
where W is the width of the phase field and H W 
is some macroscopic length scale. KR derived a kinetic 
boundary condition which relates the effective tempera- 
ture T(0) and the growth velocity V by the kinetic co- 
efficient Akr- [Tm - T{Q)]L/TIj = AkrV. Using the 
asymptotic matching procedure they obtained that the 
kinetic coefficient has the following structure: 

AKR^Tf^((3o-a—), (21) 

where /3o = 1/Vb > is the KR kinetic coefficient in the 
sharp interface limit (VF — 0) and a is a positive numer- 
ical factor of the order of unity which depends on some 
tiny details of the specific phase field model. The second 
negative term is due to the finite thickness W of the in- 
terface and the described linear extrapolation procedure. 
We also note that in this description the other Onsager 
coefficients vanish, B = C = in both sharp and and thin 
interface limits. KR checked that for the steady state 
ID growth, the analytical prediction, Eq. ([5]) with the 
obtained value of Akr and ;B = C = 0, is in good agree- 
ment with direct numerical simulations of the phase field 
model. However, there is a subtle physical point concern- 
ing the interpretation of A^ which may become negative 
with some choice of phase field model parameters. As 
correctly mentioned by KR, this conclusion may appear 
at first sight thcrmodynamically inconsistent. However, 
as it has been already mentioned, the temperature T(0) 
is not a real temperature inside of the interface and de- 
viates strongly from the temperature obtained by phase 
field simulation, which is below r(0). 

Let us discuss this nontrivial point in more details. We 
can imagine an extended interface with the thickness 25 
with two boundaries located at a; = ±5. We emphasize 



that this length scale 5 is different from the phase field 
interface width W and for the moment arbitrary, still 
being much smaller than any relevant macroscopic length 
scales. We can easily derive the corresponding matrix 
of Onsager coefficients using the values of T and /i at 
the two boundaries of the extended interface as Ti = 
T(0)"Gi(5 and T2 = T{Q) + G25, and /ii(Ti) and ^i2{T2). 
Using Eqs. ([S])-® we express the temperature gradients 
Gi in terms of Je and V, and using Eqs. ©-(H]) we finally 
find the renormalized values of the Onsager coefficients 

A{5) = Akr + C{5)TUsI + 5|)/2, (22) 
B{6)^-C{5)Tm{Si+S2)/2, (23) 
C{5)Tl, = 26/X, (24) 

where, we have assumed that Ai = A2 = A as in Q. Few 
remarks are in order. 

i) The steady state result, Eq. ([5]), is invariant with re- 
spect to this renormalization of the Onsager coefficients, 
i.e. independent of S. It means that this i5-family of 
Onsager matrixes is in good agreement with numerical 
simulations of the phase field model as well as the origi- 
nal KR case, (5 = 0. 

ii) With the choice 6 > 2aW the matrix of On- 
sager coefficients becomes positive-definite, A,C > and 
AC > B^, for arbitrary parameters of the phase field 
model. This result has a clear physical meaning. For 
6 ^ W we discuss only the range of |x| where the used 
linear extrapolation of the temperature field is in agree- 
ment with the temperature field obtained by the phase 
field, while ior S the temperature at the boundaries 
strongly deviates from the phase field description, which 
is fully thcrmodynamically consistent by itself. In other 
words, ior 6 ^ W the obtained matrix of kinetic coeffi- 
cients does describe real physical dissipation in the region 
S, while for S ^ W this "effective" matrix does not de- 
scribe the real physical dissipation, but still leads to the 
correct expression for the steady state growth velocity. 

This possible renormalization with 6, much smaller 
than any macroscopic length scale H, is not specific only 
to the phase field models and represents a small effect of 
the order of <C 1. It has the same structure as the 
"negative" phase field effects W/H. The ideology of any 
macroscopic description relies on this small parameter as 
an expansion parameter of the theory. These corrections 
should be irrelevant in the general case of the diffusional 
transformation where the bulk dissipation plays the ma- 
jor role (for example, in the case of dendritic growth at 
small undercooling). We have seen, however, that in the 
specific problem of steady state ID front propagation this 
small term (proportional to W) is responsible for the sign 
of the slope in the phase field model description. This 
happens because the bulk dissipation (being still much 
larger than the interfacial dissipation) just bring us to 
the vicinity of point At = 1 and does not contribute 
to the slope. In this case the growth velocity is entirely 
controlled by the interface kinetics. We note that the 
interpretation of the nontrivial behavior in the vicinity 
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of A = 1 due to sufficiently negative values of the phe- 
nomenological cross coefficient B does not assume any 
specific model of the interface. At the same time, the ex- 
planation suggested by the phase field modeling explicitly 
takes inhomogeneities of the temperature and concentra- 
tion fields, on the scale of finite interface thickness, into 
account. 

In other words, there is no doubt about thermody- 
namic consistency of the phase field model for arbitrary 
values of the parameter Vg W/ Dt apart from the obvious 
restrictions, Vg > and Dt > 0. However, the interpre- 
tation of the thin interface limit and its relation to the 
matrix of dissipative Onsager coefficients should be taken 
with care. We illuminate this warning by the following 
additional example. Let us assume that initially the two- 
phase system is at some temperature T slightly below the 
melting temperature Tji/. This system evolves towards 
equilibrium with a solidification velocity V that decays 
a,sV^ at large time t. This behavior would be ob- 

served in direct phase field simulations for arbitrary pa- 
rameters of the model independent of the sign of effective 
kinetic coefficient, Eq. (PT|). A slightly different but close 
in spirit non-stationary evolution has been discussed in 
Q confirming this behavior. However, if one solved this 
problem not by a direct phase field simulation but by 
solving the free boundary problem with effective bound- 
ary conditions described by the the matrix of kinetic coef- 
ficients, A = Akr and B=C=0 (the thin interface limit 
of the result would be very different if Akr < 0. 
The system would melt, instead of of being solidified, ex- 
hibiting strong instabilities and would never reach the 
described physical attractor. On the other hand, if one 
solved the same problem using the renormalized positive- 
definite matrix of Onsager coefficients, Eqs. (|22ll24p . the 
result would be basically the same as in direct phase field 
simulations and physically relevant. Therefore, we con- 
clude that the interpretation of the thin interface limit 
of 0] as the correspondence between the phase field de- 
scription and the classical macroscopic approach is in- 
correct for the wide class of non-stationary problems if 



Akr < 0. However, the renormalized positive-definite 
matrix of Onsager coefficients leads to such a correspon- 
dence in the macroscopic limit for arbitrary Akr- 

Finally, we would like to address one more point. The 
phase field model of Q contains less independent param- 
eters to describe the kinetic properties of the interface 
(only Aor /3q) than is allowed by the general phenomenol- 
ogy {A, C). While an independent parameter C can be 
introduced in a slightly modified version of the phase field 
model, the introduction of the independent cross coeffi- 
cient B is a serious problem. As pointed out in [l^, ac- 
cording to Curie's principle [l6j . there can be no kinetic 
coupling between the scalar non-conserved phase field or- 
der parameter and vectorial diffusional fiuxcs of the 
conserved quantities energy and/or concentration. Thus, 
one should not expect an independent cross coefficient 
B to appear in the effective boundary conditions, Eqs. 
pi4p and Eqs. (jlOllip . However, in the general case 
of the phenomenological macroscopic description, we do 
not doubt the existence of such a kinetic coupling at the 
interface between the normal growth velocity and normal 
diffusional fluxes through the interface. It is conceivable 
that this coupling can be introduced in modified versions 
of the phase field model where V(/)/|V0|, the unit vector 
normal to the interface, can be used to produce the cor- 
responding vectorial quantities. This issue may also be 
relevant to the anti-trapping current introduced in some 
non- variational versions of the phase field model [13, [i3| 
for different purposes. The anti-trapping current intro- 
duces a new kinetic coefficient and uses the unit vector 
normal to the interface. To use this idea for the descrip- 
tion of the cross effect of the interface kinetics in phase 
field models, one should carefully consider the necessary 
symmetry which is obligatory for this cross effect. A 
more detailed discussion of this question is far beyond 
the scope of this paper. 
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